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Abstract 


Lumped paramter LMFBR dynamical systems based on effective lifetime and prompt- 
jump approximations, and one and six groups of delayed neutrons are analyzed for 
transcritical and Hopf bifurcations. Besides the core, a simple model of secondary heat 
exchanger, as well as suitable models for control reactivity, are also included. Both the 
proportional and the differential controllers are considered. 

Data for analysis is derived for physical parameter values available for a large 
number of LMFBRs. The analysis shows that none of the LMFBRs considered is likely 
to experience Hopf bifurcation. It is also fotmd that the transcritical bifurcation occurs 
only when higher-order feedback effects are present. But even then the bifurcation occurs 
for a positive value of the coolant temperature reactivity feedback. These conclusions 
remain true even when we combine the dimensionless parameters in the worst maimer 
within certain stipulated ranges. 

Inclusion of control reactivity is shown to make bifurcations in the lumped parameter 
LMFBR models further unlikely for all negative gains. 
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Nomenclature 


A list of symbols is given below. Some symbols which are specific to a section, and whose 
meaning may be different in different sections, are not included. Standard mathematical symbols 
also do not appear in the list. 

Roman 



dimensionless 

H 


i'^co “ dimensionless 
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H 


cCpiTjg - 7],)^ / P, dimensionless 
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normal form parameter a (s) at s = 0 

b 

, a nondimensionless parameter 

H 

b(0) 
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second order coolant temperature coefficient of reactivity, K’ 
Up second order fuel temperature coefficient of reactivity, K‘ 

J3 one group delayed neutron firaction 

Ti a set of dimensionless parameters 

X one-group decay constant for delayed neutron precursors, s‘' 
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Pfb feedback reactivity 
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^ flux of thermal neutrons, m' s' 

0, 0(t) fundamental or state transition matrix 


XI 



Miscellaneous 

$ p/p 

Subscripts 

R a reference value 

0 a steady' state value 

Superscripts 

* 


a critical value 



Chapter 1 


Introduction 

1 .1 Motivation 

Mathematical formulations of physical problems often result in differential equations that 
are nonlinear. Nonlinear systems are used to describe a great variety of phenomena, in the 
physical and engineering sciences as well as in the life and social sciences. Nuclear 
reactors are no exceptions and are being increasingly analyzed as nonlinear systems. In 
particular certain concepts in bifurcation theory, chaos and fractals are likely to become 
more and more relevant to the safe operation of nuclear reactors. Bifurcations occur 
where the solutions of the nonlinear systems change in their qualitative character, as the 
parameters are changed. The theory of bifurcation, therefore, concerns all nonlinear 
systems and hence has a great variety of applications. 

It is well known that even with linear reactivity feedback effects, the dynamic models 
of the nuclear reactors are nonlinear. Nonlinearities in the reactor dynamical systems arise 
in general, due to strong interactions between the neutronics, thermal hydraulics and 
fission product poisoning. (Fission product poisoning is of less importance in LMFBRs 


which are the subject of the present study.) Under conditions involving large departures 
from the intended steady state, an intrinsically nonlinear system such as a nuclear reactor, 
may exhibit a response that has no counterpart in the linear world. 

No general method is available for solving nonlinear differential equations and. 
unlike linear systems, few nonlinear systems possess closed form solutions. The 
emergence of computers as new tools for theoretical physics and engineering smdies have 
favored the development of nonlinear dynamics. Numerical simulations have made the 
richness of nonlinear models accessible to the intuition of the physicists and engineers. 

At the present time it appears, however, that the techniques of nonlinear systems 
research have not been as fully and widely applied in the nuclear energy field as in other 
branches of science and engineering. The present work tries to bring out some aspects of 
nonlinear dynamics (local transcritical and Hopf bifurcations) in some LMFBR models. 


1.2 Objectives 

The overall objective of the present work is to study the phenomenological lumped 
parameter models of LMFBRs in order to understand the role of reactor core 
nonlinearities arising due to fuel and coolant temperature reactivity feedbacks. A simple 
model of the secondary heat exchanger (HE) and control reactivity feedback is also 
included. Emphasis will be on identifying the parameter ranges over which the local 
bifurcation phenomenon may take place. 


1.3 Review of Literature 

The literature in the area of LMFBR dynamical systems is reviewed here. A tvv^o 
temperature feedback model for the PWR is discussed by Manmohan in his doctoral 
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thesis (1996). Based on the above work a two temperature feedback model for the 
LMFBR is presented in this work. 

The mathematical and physical importance of the nonlinear form of reactor kinetics 
equations, when the effect of temperature is included, was pointed out by Chemick 
(1951). Robinson (1954, 1955) extended the treatment of Chemick (1951) in regard to 
concept of reactor stability and included the effect of delayed neutrons. The concepts of 
orbital, secular and asymptotic stability are examined in the phase plane representation 
and techniques are pointed out whereby the effect of the delayed neutrons on reactor 
kinetics may be taken into consideration. Hsu (1968), Kastenberg (1968), and Hsu and 
Sha (1969) considered the stability problem for spatially dependent nonlinear reactor 
systems and pointed out that in a practical reactor the delayed neutrons must be included 
and more than one feedback should be considered. Enginol (1985) has analyzed 
asymptotic stability of nuclear reactors with arbitrary nonlinear feedback. Adebiyi and 
Harms(1989, 1990) examine some aspects of linear and nonlinear nuclear reactor kinetics 
by general topological methods. Doming(1989) provides a general discussion of strange 
atrractors in the context of nuclear reactors. Yang and Cho (1992) present expansion 
methods for finding nonlinear stability domains of nuclear reactor models. 

A model for simulating the transients in LMFBR is discussed by Agrawal and 
Yang(1977). Lewins (1987) discusses various methods for the control of nuclear reactors. 
The models for the secondary heat exchanger are also discussed. Edwards and Lee (1990) 
have described a new control concept, state feedback assisted classical control for 
nuclear reactors and power plants. An improved reactivity table model for LMFBR 
dynamics is described by Stacey (1972). Some of the physical, technical, safety, and 
material problems of fast reactors are presented by Karl Wirtz (1973). 

. Hummel and Okrent (1970) have discussed the reactivity coefficients and their effects 
in LMFBRs. Data and methods for calculating reactivity coefficients are introduced. The 
significance of reactivity coefficients for the operation and safety of LMFBRs is well 
discussed. A detailed discussion of transcritical and Hopf Bifurcations is provided in 
Guckenheimer and Holmes (1983). 
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1.4 Outline of the present work 


A brief description of LMFBRs and reactivity effects of temperature is presented in 
section 1 .5 later in this chapter. Chapter 2 consists of formulations of lumped parameter 
models for the LMFBR dynamical systems for the core, secondary heat exchanger and 
control systems based on one or Six groups of delayed neutrons. Models based on 
effective lifetime and prompt-jump approximations are also presented. All the governing 
equations are presented in a dimensionless form. 

In Chapter 3, the phenomenon of transcritical and Hopf bifurcation in the LMFBR 
dynamical systems is investigated. Center manifold is calculated for the transcritical and 
Hopf bifurcations . Various numerical methods involved in calculating the bifurcation 
points in different models are also discussed. 

The behavior of LMFBR dynamical systems in the close vicinity of the bifurcation 
point is studied by numerical simulations. The results of these are presented in Chapter 4. 
A summary of work done and the possible extension of the present work, are also 
presented in Chapter 4. 

1.5 A Brief Description Of LMFBRs. 

If the main source of fission is the capture of fast neutrons directly by the fuel without the 
neutrons having suffered appreciable energy losses, the system is called a fast neutronic 
system. The fuel consists of two parts, a core, (highly) enriched in U and a 
U“*/Th232 breeder blanket. The reactor is sodium cooled and operates at temperatures 
which do not approach boiling. Such reactors are called Liquid Metal Fast Breeder 
Reactors. 

The use of liquid sodium, Na-23, as coolant ensures that there is little neutron moderation 
in the fast reactor. The element sodium melts at 37 IK, boils at 1156K, and has excellent 
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heat transfer properties. With such a high melting point, pipes containing sodium must be 
heated electrically and thermally insulated to prevent freezing. The coolant becomes 
radioactive by neutron absorption, producing the 15h halflife Na-24. Great care must be 
taken to prevent contact between sodium and water or air, which would result in a serious 
fire, accompanied by the spread of radioactivity. To avoid such an event, an intermediate 
heat exchanger is employed, in which heat is transferred from radioactive sodium to 
nonradioactive sodium. 

Some of the major differences between thermal reactors tind fast breeder reactors, from a 
dynamical viewpoint are: 

(1) The prompt lifetime of the fast reactors is a few orders of magnitude smaller than that 
of thermal reactors, being of the order of 10"’ s. 

(2) The temperatvire coefficients of reactivity are generally smaller than those in the 
thermal reactors and there is no fast acting, large dominant temperature coefficient which 
overrides all the others. 

(3) The delayed neutron fraction p for fast fission of is about fifty percent smaller 
than that for the thermal fission of U . 

(4) The total excess reactivity in a fast reactor is small. The excess reactivity is small 
because fission-product poisoning is small, and the amount of reactivity needed for 
temperature override is small. 
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Figure 1: A schematic diagram of a two-temperature LMFBR 
model with secondary heat exchanger 
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1.6 Reactivity effects in LMFBRs 


In a LMFBR three effects on reactivity require special attention which are: 

i) Sodium void effect, 

ii) Doppler effect, and 

iii) other temperature effects. 

Although the reactivity effect due to fission product poisoning of xenon and samarium 
plays an important role in thermal reactors, we can safely neglect their effect in a fast 
reactor such as an LMFBR, as the absorption cross-section of xenon and samarium is 
very low for fast neutrons. We will now briefly explain the reactivity effects due to 
sodium voiding , Doppler broadening of fission and capture resonances, and temperature. 


Sodium void effect occurs due to loss of some or all sodium fi:om a fast reactor. 
Reactivity temperature coefficients in a fast reactor can be both prompt and delayed. The 
prompt coefficient arises mainly from the Doppler effect in the fuel, but in some cases 
neutron leakage can be important. The delayed coefficient results firom a number of 
factors associated with the expansion of the liquid sodium coolant when the temperature 
is increased. 

The reactivity effect resulting fi'om Doppler broadening of fission and capture 
resonances arises firom a varying competition among fission, capture, and leakage 
processes. The Doppler effect (by Doppler effect we mean the change in multiplication 
factor associated with an arbitrary change in fuel temperature) in large, fast power 
reactors occurs almost entirely below 25 kev, because cross section variations with 
temperature become very small at higher energies . Since leakage in such reactors is 
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relatively small at low energies, their Doppler effect can be regarded as resulting from a 
varying competition between fission and capture. This phenomenon causes a difference in 
fission-source neutrons produced per neutron absorption at low energies. 

The change in this ratio (neutrons produced per neutron absorbed) for some energ>' 
below which practically all the Doppler effect occurs (perhaps 25 kev) multiplied by the 
fraction of total fission-source neutrons slowing down below that energy is approximately 
equal to the Doppler effect. The importance of the Doppler effect lies in its provision, for 
reactors with a high concentration of fertile material, of a prompt, negative reactivity 
effect from fuel heating. 
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Chapter 2 


Lumped Parameter LMFBR 
Dynamical Systems 

2.1 Neutron Kinetics 


The neutron kinetics equations with one group of delayed neutrons can be written as 
(Lewins 1978, Hetrick 1965) : 


dN 

dt 


^-^N+;iC, and 
A 


( 2 . 1 ) 


^ = Zn-AC, (2.2) 

dt A 


where 

N = N(t) is the number of neutrons or average neutron density in the reactor. 




C = C(t) is the number of delayed neutron precursors or average delayed neutron 
precursor density, 

t = time (s), 

A = neutron reproduction or generation time (s), 

P = delayed neutron fraction (dimensionless), 

X = decay constant for delayed neutron precursors (s"‘ ), and 
p = reactivity or fractional excess multiplication constant (dimensionless). 

2.2 Power Balances 

The two temperature feedback model is based on two power balances: one for the 
energy contained in the fuel elements of the reactor core, and another for the energy 
stored in the coolant within the core volume. Although a small fraction of energy released 
in the fission process is directly deposited in the coolant, this is ignored. 


The power balance for the fuel elements Is given by 


dT. 


/ _ 


dt 


P-H,(T, -TJ, 


( 2 . 3 ) 


where 
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T j = (t) is the average fuel temperature(K), 

T -^T 

T ^ = T ^ (t) = is the mean coolant temperature(K), 

P = P(t) is the reactor pov.'er(W), 

C ^ = heat capacity of the fuel elements in the reactor core(J K'^), 
and 

Hy = fuel to coolant heat transfer constant(W K’^). 


The power balance for the coolant in the reactor core is given by 

=H,fr/-TJ-2mc.(T, -T,). (2.4) 

where 

C ,. = heat capacity of the coolant in the reactor core(J K'*), 

T, = coolant inlet temperature(K), 

/w (. = mass flow rate of coolant through the core (kg s"'), 

Cg = specific heat of the coolant (J kg'' K''). 


The power balance for the secondary heat exchangens given by 


c, ^ = H,(T, - r,) -2m, cXT. - z,) 


(2.5) 


where 


C j = heat capacity for the secondary side(J K' ), 

= secondary side inlet temperature (K), 

T +T- 

Tj = — — is the average temperature corresponding to the secondary side (K), 
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= Secondary side flow rate(kg s"') 

= primary to secondary side heat transfer constant (WK'*) 

2.3 Steady State Relations 


The steady state operating values of the state variables in Eqs. 2. 1-2.5 will be denoted by 
N„ , C„ , T Tso and respectively, and the steady state value of the reactivity is 
taken as . By equating Eqs 2. 1-2.5 to zero, we obtain 


Po = 0, 


Cq 



o 



2mc^ 



) 




( 2 . 6 ) 


Tjo is the steady state value of the reactor inlet temperature, Tj. It is assumed to be 
constant if secondary heat exchanger is not included in the model. Otherwise, Tj;, the 
inlet temperature to the secondary heat exchanger is considered constant. 


2.4 Linear reactivity Feedback 


The reactivity p appearing in Eq. 2.1 can, in.general, be written as 
P Pa + A(t)+Pc+Pyi 


(2.7) 
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As seen above, if there is no external source of neutrons in the reactor, as is the case 
considered throughout the present work, = 0. For the autonomous systems we treat in 
this work, ~ 0- The control reactivity p^ is considered separately in section 2.10 
Further, the dynamical systems developed in this section are based on linear reactivit>‘ 
feedback which for the two temperature LMFBR model can be written as 

P = Pfi = (Xf(J f ~ Tj.J+q;,(T^ - T,.„), (2.8) 

where 

= fuel temperature coefficient of reactivity(K ), and 

= coolant temperature coefficient of reactivity(K"‘). 

Higher order feedback terms are considered later in section 2.8. 

It is the substitution of Eq. 2.8 in Eq. 2.1 which is one source of nonlinearity in the 
fission reactor dynamical systems. It leads to the product terms of N(t) or P(t) with both 
T j (t) and Tc(t). Instead of presenting the resulting equations, we first choose suitable 
nondimensional forms of the dependent state variables and the independent time variable. 

2.5 Dimensionless Equations 

In order to obtain the previous equations in a dimensionless form we choose the 
following translation and scaling of the state variables and time: 


P-P„_N-N„ _ C-Q 

_ . _ 5 ^ 


III 


1 



1 

Fh" 

II 

Tp-T/ 

c 


—t , and 

= 

r 

1 

A 


T -T 

* so 


(2.9) 
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where Tj = Tjo if the secondary heat exchanger is excluded from the model, and Tj = Tjo, 
if the secondary heat exchanger is included. 


The translations of the state variables in the above equations shift the operating or 
equilibrium point to the origin. This is essential for die analysis presented in later 
chapters. The normalizations used could be chosen differently, but the above choice puts 
the dynamical system in the simplest form. Substitution of Eqs. 2.8 and 2.9 into Eqs. 2.1- 
2.5 and using Eqs. 2.6 where needed leads to the following dimensionless system: 


dp 

dr 

= -p+c+a^ 3^+a, 3,+a^ p3„ 

(2.10) 

d d 

dr 

u 

1 

II 

(2.11) 

d3j- 

dr 

= p(l-0)p-p3^+p53,. 

(2.12) 


If the secondary heat exchanger is excluded from the model, we have for the 
dimensionless coolant temperature 3 ^ : 


dr ^ ‘ 


(2.13) 


If the secondary heat exchanger is included in the model then we have the following 
equations for 3 ^ and 3^ : 


dr 


= E.{ 

d ^ 




(2.13a) 


d3 

dr 


= r(3,- (1 + 0)3,), 


(2.14) 


where 




Sc = 




14 



b 


P 


AA 


P ’ 


APq 

pC,{Tj,-T,,y 


G = 



c 





r 


^0 

PCAT„-T,)’ 


T -T 

© = _£2 *2. (2.15) 

rp rp N ^ 


It can been seen that the original system of two temperature feedback model has been 
transformed to a dimensionless system with six or eight nondimensional parameters, 
depending on whether the secondary heat exchanger is excluded or included in the model. 
The parameters a^- and a^ can been interpreted as the dimensionless measures of 
temperature coefficients of reactivity. The parameter p may be understood to measure the 
steady state power Po of the reactor. The other parameters are viewed as simple 
dimensionless combinations of input quantities involving times, temperatures and heat 
capacities. 

Equations 2.10-2.14 can been written in the form of a state equation: 

X = V(X,Ti) 

= U(Ti)X + W(X,q) (2.16) 

where the overdot denotes the derivative with respect to the dimensionless time r , X is 
the state vector, r\ is the set of dimensionless parameters (Eqs. 2.15), and V(X,ri) is the 
vector field whose linear part is U(ti)X and the non linear part is W(X,ti), i.e., 


1 


15 





' P' 


+a^X^ 



a 


0 


= 


, W(X,T1) = 

0 

^4 


5, 


0 

.^ 5 . 


.3.V. 


0 



-1 

1 



0 


b 

-b 

0 

0 

0 

U(T,) = D,V(X„,ti) = 

pO--e) 

-P 

pe 

0 

0 


0 

0 

d 

-d 

0 


0 

0 

0 

r 

-r(l + 0) 

r] = {a^ b c p 

» r © } 






P 


-1 

1 



0 


' P' 

c 


b 

-b 

0 

0 

0 


c 

2/ 

== 

p{\-e) 

-P 

pe 

0 

0 



3, 


0 

0 

d 

-d 

0 





0 

0 

0 

r 

-r(l + 0) 

1 



2.6 Six Group of Delayed Neutron 
Precursors 

Six groups appear to be optimum for an accurate representation of delayed neutron 
precursors. For the calculation of the critical value of a parameter at which a bifurcation 
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may take place, it is therefore desirable, to use six groups of delayed neutrons. This is 
achieved by replacing Eq. 2.1 1 by six differential equations governing the population of 
different delayed neutron precursors: 


^=b,( p - c,), 


i=l,...,6, 


where 


<=; - 


c:, = population of the precursor group i, 
cZiQ = steady state population of the precursor group i, 
2., A 

' X 


Pi = delayed neutron fraction for the ith precursor group, 
decay constant of the ith precursor group, 

and other variables or constants are as defined earlier in this chapter. The term e in Eqs 
2.10 is replaced by: 


/a=l 

where a, is the relative yield for the precursor group i. For the purpose of 

parameter variations in six group models, instead of treating each bj as a variable 
parameter, we write 


, j , 1 

b: = — 0 , and a = 

' \X) 

jLa ^ 
/«i 
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where b has been previously defined (Eqs 2.15) and is treated as a possibly variable 
parameter as in one-group models. 


2.7 Simpler Dynamical Models 


In the dynamical models presented in the previous section, the delayed neutrons have 
been treated using a one or six precursor groups. One precursor group is often convenient 
and satisfactory. Further simplification is afforded by the effective lifetime and the 
prompt jump approximations. Below these models are briefly described and the 
governing equations are presented in the dimensionless form as used in the present work. 


Effective Lifetime Model 


1 A 


The effective lifetime A of neutrons may be defined as (Hetrick, 1971 ; Lewins, 1978): 

(2.18) 

If, in the preceding equations, A is replaced by A^, and the state variables as well as the 
equations corresponding to the delayed neutron precursor are dropped, we obtain the 
lumped parameter dynamic equations based on the effective lifetime of neutrons. This 
model is sometimes easier to analyze but is applicable only when the reactivity and the 
rate of change reactivity are small : 


dt 


«\A « 


(2.19) 


or, in terms of the dimensionless variables used in this work, 

u$i 


$ 1 «1, and 


dr 


«bl$| 




where $=~ and b is defined earlier in this chapter. 
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Prompt-Jump Approximation 


The prompt-jump approximation or, as it is sometimes called, the zero-lifetime 
approximation, is described in detail by Hetrick (1971). Formally, an expansion of the 
neutron population or density in powers of the small parameter A is attempted, the 
procedure being analogous to the method of singular perturbations. If this approximation 
is used, the point-reactor kinetics Eqs. 2. 1 and 2.2 are replaced by the following single 
equation : 


dt 


Xp + 


dt 


P-p 


N 


( 2 . 20 ) 


In Eq. 2.20 the reactivity p must be substituted through Eq. 2.8, in terms of the other 
variables. This makes Eq. 2.20 for the neutron population nonlinear, as is the case in 
other models. The criterion for the validity of prompt-jump approximation is (Hetrick, 
1971): 



dN 

«P-P 

N 

dt 

A 


( 2 . 21 ) 


It should be noted that the criterion (2.21) can only be satisfied it p < p , i.e., $ < 1. In 
terms of the dimensionless variables defined earlier, the prompt-jump approximation can 
be written as 



(l-^p)($-^6S) 


( 2 . 22 ) 


which is applicable if the following condition, equivalent to (Eq. 2.21), is satisfied: 
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« 1 . 


(2.23) 


P 

(l + p)(l-$) 

where the overdot denotes derivatives with respect to the dimensionless time r defined 
earlier. 


2.8 Higher-Order Feedback 


Higher-order feedback terms do not play a significant role in Hopf bifurcation 
(Manmohan, 1996). However, as will be seen later in this work, they play a crucial role 
in the transcritical bifurcation. We describe below the higher-order feedback reactivity as 


follows: 

P = Pfb = Pft + Pfb^^^ + > 

(2.24) 

with 





(2.25) 


Ptt® = -Tj,r+ ac(J, -TJ\ 

(2.26) 


where ap and denote the second order reactivity feedback coefficients for the fuel and 
the coolant , respectively, and the other symbols are as used earlier in this chapter. It is to 
be noted that the units of ap and ac are K'^, while those of the linear feedback 
coefficients, a jr and are K"'. Introducing Eqs. 2.24-2.26 in Eqs 2.7 and 2.1, and 
making use of the nondimensional variables and parameters already defined, we can write 
the first equation in reactor kinetics as follows: 


dr 


= -^+c:+a^^j- 




(2.27) 
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where 


Op =ap{Tj,-TjY I 
Qp — OCp{T^g ^ Pj 

and all other parameters are as defined earlier. 


2.9 Control Reactivity Feedback 


In the present work we consider a simple control model based on proportional and 
differential controllers. The control may be based on the fi-actional change in reactor 
power and / or a suitably defined fi'actional change in reactor coolant temperature. Thus 
the control reactivity can be written as 


where 

= 

K = 




( 2 . 28 ) 


(negative) constant (gain) proportional to the fractional change in power 
represented by the dimensionless power p 

(negative) constant (gain) proportional to the firactional change in coolant 
temperature represented by the dimensionless coolant temperature 3^. , and 
time constants associated with the differential controller based on the rate of 
change of power or coolant temperature, respectively. 


IncltKion of control reactivity changes only the first equation in the original system 
presented earlier. If the differential controller is included then the resulting equation is 
suitably recast by transferring the derivative terms on the LHS, or substituting their 
values in terms of state variable from other equations. ‘‘ 


The resulting equation governing the dimensionless power is given be] 
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1 


\-kr 

P P 




p = (-l + -fcp)p+c+(fl^ ■¥kj^d):5y +(a, -k^T^d):5^ 


where 




A 


, and T = 


A 

A 


Here only the linear terms on the RHS have been retained. As we will see later, analysis 
based on these linear terms show that bifurcation (static or dynamic) takes place for very 
large (negative) values of or , which are unlikely to occur in practice. As a result 
the higher-order terms on the RHS are not required in the present work. 

2.10 Input Data 


The set of input values that we need to know for each reactor model to carry out the 
analysis and computations have been arrived at in the preceding sections of this chapter. 
The input data are evaluated from the reactor data available from Power Reactors, 
(1983)and World Nuclear Industry Handbook (1994). Data of eleven different LMFBRs 
is given in Table 1. For reactors of different designs, or under different operating 
conditions, the input data can be widely different. 



Table 1 : Input data ' for LMFBRs 


Reactor 

Po 

Cf 

Cc 

Tfo 

Too 

T- 

MWt 

MJK'' 

MJK'* 

K 

K 

K 

SUPER-PHENIX 

3000 

5.355 

8.21 

1583 

743.0 

668 

BELOYARSK-3 

1470 

2.72 

1.538 

1628 

736.5 

650 

BELOYARSK-4 

2100 

2.72 

1.538 

1628 

736.5 

650 

DFR 

60 

0.057 

0.09 

1283 

754.0 

673 

EBR-2 

62.5 

0.078 

0.07 

868 

694.5 

644 

JOYO 

100 

0.08 

0.186 

1848 

708.0 

643 

KARLSRUHE 

58 

0.232 

0.226 

1865.5 

715.5 

633 

MANGYSHLASKI 

1000 

0.39 

1.497 

1373 

628.0 

553 

MONJU 

714 

1.888 

1.635 

1785.5 

736.0 

670 

PHENIX 

563 

1.376 

0.773 

1773 

753.0 

673 

FBTR 

40 

0.047 

0.03 

1473 

718.5 

653 


t The meaning of each symbol is given in the Nomenclature and explained in the text. 


A sample calculation for the values of Cf and C,. is given in Appendix A. Table 2 gives 
the reference values for one group reactor-kinetics parameters as well as the 
corresponding dimensionless parameters. Table 3 gives the six group of delayed neutron 
data. The dimensionless parameters are given in Table 4 are calculated based on data 
given in Table 1. The value of the parameter r in Table 4 is based on, 

Cj = 0.02 Full-power-second (FPS). 
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Table 2: Input data and dimensionless parameters for one-group reactor-kinetics 


1 




Input Values 

Dimensionless Parameters 

Symbol 

Value 

Symbol 

Value 

P 

0.0035 

b 

2x lO"^ 

X 

0.07 s'' 

Sf 

-1.045 

A 

10’’ s'' 

ac 

-0.128 

Uf 

-4x 10'^ K'' 

t 

r 

0.000029 


-6x 10*^ K'' 




f The meaning of each symbol is given in the Nomenclature and explained in the 
text. 

I This value of r is calculated on the basis of : Cj = 0.02 FPS, 

T3o = T,o-50,T,i = Tio-50. 


Table 3: Input data for six groups of delayed neutrons 


Group 

i 

Relative Yield 

at 

Decay Constant 

X, (s-‘) 

1 

0.038 

0.0127 

2 

0.213 

0.0317 

3 

0.188 

0.115 

4 

0.407 

0.311 

5 

0.128 

1.40 

6 

0.026 

3.87 
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Table 4.: Dimensionless parameters for LMFBRs based on data given in table 1 


Reactor 

P 

C 

e 

1 

0 

SUPER-PHENIX 

0.000019 

0.652 

0.081 

0.667 

BELOYARSK-3 

0.000017 

1.768 

0.088 

0.578 

BELOYARSK-4 

0.000025 

1.768 

0.088 

0.578 

DPR 

0.000057 

0.633 

0.132 

0.617 

EBR-2 

0.000132 

1.114 

0.225 

0.990 

JOYO 

0.000031 

0.43 

0.053 

0.769 

KARLSRUHE 

0.000006 

1.026 

0.066 

0.606 

MANGYSHLASKI 

0.000098 

0.26 

0.091 

0.667 

MONJU 

0.000010 

1.154 

0.059 

0.757 

PHENIX 

0.000011 

1.78 

0.072 

0.625 

FBTR 

0.000032 

1.566 

0.079 

0.763 


t These values are calculated on the basis of : Cj = 0.02FPS, 

T = T -50 T • = T- -50 

•*■30 •*■00 ■*■31 •‘•10 

t The meaning of each symbol is given in the Nomenclature and explained in the 
text. 
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Chapter 3 


Analytical and Computational 
Techniques 


Several analytical and computational techniques are needed to locate the bifurcation 
points and to investigate the behavior of the dynamical system in the vicinity of these 
points. These have been discussed in detail by Wiggins (1989) , and Parker and Chua 
(1989 ) . The techniques required for the analysis of Hopf bifurcation are also developed 
by Manmohan (1996). A bifurcation describes a qualitative change in the dynamics when 
the parameters in the system are varied. The values of the parameters where this change 
takes place are called bifurcation values or critical values. Knowledge of bifurcation is 
absolutely essential for the complete understanding of the system. Below we summarize 
the methods used in the present work for the study of both the static and Hopf 
bifurcations in LMFBR systems. 

3.1 Normal Forms at Bifurcation Points 

A static bifurcation occurs when an eigenvalue crosses the imaginary axis through the 
origin. On the other hand if a pair of complex conjugate eigenvalues cross the imaginary 



axis transversally a Hopf bifurcation occurs This is illustrated in Figures 2a and 2b. In 
both the cases three different situations may arise. For Hopf bifurcation this is illustrated 
in Figures 3a-3c, and for the static bifurcation this is illustrated in Figures 4a-4c. 
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Figure 2 : A schematic diagram of (a) static bifurcation (b) hopf 
bifurcation 
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Figure 3: Schematic diagrams of (a) saddle-node (b) pitchfork 
and (c) transcritical bifurcations 



Parameter e 



Figure 4: (a) A schematic diagram of a normal (Supercritical) 
Hopf bifurcation 



Figure 4 : (b) A schematic diagram of a subcritical (inverse) 
Hopf bifurcation 




X 



Figure 4 : (c) A schematic diagram of Coexistence of unstable 
and stable solutions in a subcritical Hopf bifurcation 


I-' 

k\. 
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Poincare Normal Form 


A systematic way to ascertain the effect of the nonlinear terms is to reduce the dynamical 
system to a two dimensional center manifold, which is then transformed to the Poincare 
normal form. The proof that this can be done may be found in Wiggins (1989). Here, we 
state the main results related to the Poincare normal form and introduce some notations 
used in the calculations. The Poincare normal form in the Cartesian coordinates, Xi and Xi 
for Hopf bifurcation, can be written as 


Xj = a{s)x^ -co{s)x 2 + [a(s)x^ +^ 2 ^)+ P-l) 

Xj = 4u(5)x, + a(s)x2 + [i(fi')Xi + a(s)x2 ]ix■^^ + Xj^ ) + (3.2) 


where the normal form parameters a(s) and b(e) are defined for all sufficiently small s. 
The parameter a(0) , called the stability parameter, determines the stability of 
bifurcating periodic solutions (Wiggins 1989); 

( i ) The case <af(0) < 0 : The bifurcation is normal (supercritical). 

( ii ) The case a(0)> 0 : The bifurcation is inverse (subcritical). 

Remark : A simple proof of the above statements, facilitated by considering Poincare 
normal form in polar coordinates, may be found in Wiggins (1989). 

We now give below the normal forms for the three kinds of static bifurcations namely 
saddle-node, pitchfork, and transcritical bifurcation 

• The normal form for the saddle-node bifurcation is given by, 
x = e±x^ 

• The normal form for the pitchfork bifurcation is given by. 
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X = sx + x^ 

• The normal form for transcritical bifurcation is given by, 
x = sx + x^ 

In the above equations e is the variable parameter. Let p, e 77 , and p* be the value of p at 
which the bifurcation occurs. Then s is defined as 

s = tzJL 

3.2 Location of Bifurcation Points 


Several techmques have been discussed in the literature for the location of the bifurcation 
point (Kubicek and Marek, 1983). In general the algorithms used are iterative in nature 
and make use of the Newton-Raphson method to find the critical value of the parameter. 
The critical value p* of p e 77 is located by solving the equation 

a (p) =0 (3.3) 

where a (p) denotes the real part of a complex conjugate pair of eigenvalues of the 
matrix U(ri), evaluated at the steady state / equilibrium (operating) point. 

The QR Method 

The eigenvalues of the matrix U(ti) for a given set of parameter values are calculated 
using the QR algorithm. The Jacobian matrices associated with the dynamical systems 
studied in the present work are, in general, nonsymmetric. Since a nonsymmetric matrix 
may not be balanced, it is first replaced by a balanced matrix with identical eigenvalues 
(Press, et al., 1993). This reduces the sensitivity of the eigenvalues to rounding errors 
during numerical computations. The matrix is then reduced to a simpler Hessenberg form 
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and the eigenvalues are calculated by using the standard QR algorithm applicable to real 
Hessenberg matrices. 


Newton-Raphson iteration 

The Newton-Raphson algorithm is used to find the solution of Eq 3.3, Since the success 
of a locally convergent scheme depends critically on having a good first guess for the 
solution, it is desirable to use a globally convergent method which guarantees some 
progress towards the solution at each iteration. Perhaps the best strategy is to bracket the 
root and then refine it by a combination of the Newton-Raphson and the bisection 
methods. This keeps the root within the brackets while taking advantage of the rapid local 
convergence of the Newton-Raphson method. We have mostly used this kind of approach 
in this work. 

3.3 Calculation of Center Manifolds 


Center manifolds play an important role in the organization of system orbits in the phase 
space. This is due to the fact that these are invariant under the flow generated by the 
dynamical system. The center manifold plays a particularly crucial role in the study of a 
bifurcation. It isolates the complicated asymptotic behavior by locating the invariant 
manifold tangent to the center eigenspace, i.e., the subspace spanned by the eigenvectors 
corresponding to the eigenvalues on the imaginary axis. 


Hopf Bifurcation 

Calculation of center manifolds and normal forms requires considerable amount of 
algebra, [referred to as Horrendous by Wiggins (1989) ] particularly for the 2- 
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dimensional Hopf bifurcation. These algebraic aspects are investigated in detail by 
Manmohan (1996) and in the present work we used the same computer program to 
determine the stability parameter, a(0), occurring in the Poincare normal form. We first 
consider the parameter independent case and let U=U(rj’) and W(X) = W(X,ri*) to 
simplify the notation. 

In order to calculate the center manifold of the dynamical system 

X = UX + W(X) 
let X= SY 

where S is the matrix whose columns are the eigenvectors of U. Thus we have 

Y=(S'‘US)Y + S’‘W(SY) 

= AY + S‘‘W(SY) 

Equation 3.2 can be separated into tw'o parts : 

u = Aw + F(u,v) 

V = Bv + G(u,v) 

where 



0 Q)„ 


u 

A = 


, Y = 



_-Q}„ 0 _ 


V 


Then the introduction of another nonlinear transformation 

V = h(u) (3.8) 

followed by elimination of quadratic terms in the resulting equations is required to obtain 
the Poincare normal form. For details reference should be made to Wiggins(1989), 
Manmohan( 1 996). 

Static Bifurcation 

The transformations required for the calculation of center manifold/nonnal form for static 
bifurcation are similar to those for the Hopf bifurcation, but considerably simpler due to 


(3.4) 

(3.5) 

(3.6) 

(3.7) 
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one-dimensional nature of the manifold. We illustrate it using one-group LMFBR model., 
The center manifold u, (Eq. 3.7) is one-dimensional and we can write it as 


u = F( u,v) 

let Vi = hi(u), i=l,2 ,n-l. 

= bjU- 

As shown in Wiggins(1989), each hi(u) must satisfy 


dh 


du 


^F(u,h(u))-Aihi -Gi(u,h(u)) = 0, i = l,2,...n-l. 


(3.9) 


where 



‘0 . 

. . . O' 


. A, 

. 

A=. 

. 

yV.2 • • « 


_0 . 



is a diagonal matrix. 


It is straightforward to see that if only linear feedback is used solution of Eq 3.9 yields, 

b; = 0, i=l,...,n-l. 

Thus, the center manifold / normal form has no quadratic terms. If all bj’s are zero, and 
the original nonlinear terms in W(X) do not contain cubic terms (as is the case with linear 
feedback), the center manifold or normal form also do not contain any cubic terms. 

Thus with linear feedback, all the three types of static bifurcation are ruled out, and 
the behavior of the nonlinear system remains divergent as would have been the case in 
the linear system. This is in conformity with the fact that, with linear feedback, no 
additional equilibrium (fixed) point originates from the operating point, i.e., (0,0,. ..0). 

If higher-order reactivity feedback effects are present, then, as shown in Section 2.8, 
these introduce additional quadratic and cubic terms in the system (Eq. 2.27) , Due to the 
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cubic terms, an additional equilibrium solution / fixed point of interest emerges. It is easy 
to see that without the secondary heat exchanger, this fixed point is given by (using one 
group of delayed neutrons) 


=3, =- 


(X p Q.^ 


and with the secondary heat exchanger, this is given by 


P =c= 0 , 


where 


^ + 0 ^ 

=— 3,, 




1 + 0 


3 . = 


' 1 + 0 ’ 


3, = 


{-ka + A) ± SQRTiika + Af - 4kaA) 
2kA 


where 


SQRT denotes the square root. 
© 


k = 


a = 


A = 


1 + 0 ’ 

(9 + 0 

a. 

1 + 0 ^ 


af +a.. 


9+e 

1 + 09 


Of +a^ 


(3.10) 


In Eq 3.10 only one of the solution bifurcates from the origin, indicating transcritical 
bifurcation, and excluding the possibility of saddle node and pitchfork bifurcations. (The 
saddle node bifurcation is also excluded by the fact that the original equilibrium point, 
(0,0,...0), does not remain stable beyond the bifurcation point). The fact that at the 
bifurcation point, exchange of stability occurs, is brought out in the diagrams shown in 
Chapter 4. 
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As can been seen from Eq 2.27, the higher-order feedback effects introduce into the 
system not only the cubic terms, but additional quadratic terms. As a result, the quadratic 
terms in the normal form can no longer be zero and thus the normal form will correspond 
to that of transcritical bifurcation. 


3.4 Calculation of Trajectories 

Calculation of trajectories by numerical integration is perhaps the most important and 
expensive task in the simulation of a dynamical system. It is also a prerequisite for the 
application of geometric tools or mathematical constructs such as Poincare sections, 
phase portraits and the calculation of characteristic multipliers or Lyapimov exponents. 

Over the years, a large number of algorithms have been developed to approximate the 
behavior of a continuos-time system on a digital computer. A proper choice of the 
method and a careful application play a cmcial role in the computational accuracy and 
efficiency. Among the available methods, the major practical choices are the single step 
methods, the multistep predictor-corrector methods, and the Bulirsch-Stoer extrapolation 
methods. The most popular among the single step methods are the Runge-Kutta methods. 
Both explicit and implicit implementations of the methods mentioned above are available 
in literature. The choice depends on the stiffiiess of the system to be integrated. In the 
present work an implicit implementation of Bulirsch-Stoer method given by Press, etal 
(1993) is used. 
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3.5 Poincare sections and Characteristic 
Multipliers 


In order to simplify the phase space diagrams of a dynamical system of order n, one 
introduces an (n-1) dimensional surface of sections in the phase space and. instead of 
studying a complete trajectory, one monitors only the points of its intersection with this 
surface. A set of points of intersection of the trajectory with a hypersurface, as the 
trajectory crosses the surface from one side to the other, is referred to as Poincare section. 
Mapping an intersection point onto the subsequent intersection point is referred to as the 
Poincare (return) map. 

The stability of periodic solutions or closed orbits is determined by the characteristic 
(Floquet) multipliers or roots. These can be regarded as a generalization of the 
eigenvalues at an equilibrium point. According to standard Floquet theory, the stability of 
a periodic orbit is determined by the eigenvalues of the monodromy matrix ‘5(T) ( also 
referred to as Floquet transition matrix), which is simply defined as the value of the 
fundamental matrix at the time period T, with the initial value as the unit matrix I 
(Kubicek and Marek 1983). If these eigenvalues lie within the unit circle in the complex 
plane, the periodic solution is stable. It is well known that one of the eigenvalues of 
<E>(T) is always unity (Husseyin, 1986; Parker and Chua 1989). Thus a periodic orbit can 
never be asymptotically stable. The remaining eigen values of 0(T) determine the 
orbital stability and are called the characteristic multipliers. 

In the present work, we confirm the existence of limit cycles (whenever these occur) 
using Poincare sections, and calculate the characteristic multipliers using a method based 
on Poincare return maps (Manmohan 1996). The numerical results for both the 
transcritical and Hopf bifurcations are reported and discussed in the next chapter. 
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Chapter 4 


Numerical Results, Discussions, 
and Conclusions 


Numerical simulations are required to confirm the predictions of the theoretical 
considerations given in Chapter 3. Furthermore, the theory is valid only in an arbitrarily 
small neighborhood of the bifurcation points. As one moves away from this 
neighborhood, numerical calculations are often the only means of obtaining information 
about the nature of solutions. In this chapter, we present the numerical results for the 
lumped parameter LMFBR dynamical s>’stems developed in Chapter 2. 

4.1 stability at Full Power Operation 

Before proceeding to calculate the critical values of parameters for which bifurcations 
may occur and simulating the behavior subsequent to bifurcation, it is desirable to 
confirm the stability of the systems considered at their normal operating values. In reactor 
dynamical systems this is indicated by reactor period, Tp, or the e-folding time 
(Glasstone, 1994) 


where cOq is the real part of the eigenvalue closest to the imaginary axis. The eigenvalues 
are calculated for the Jacobian matrix, U (Eq. 2.17) , evaluated at the normal operating 
point. A negative reactor period indicates stability at the point under consideration. 


Table 5 shows the (negative) reactor periods for the reference LMFBR, SUPER 
PHENIX, using five dynamical models presented in Chapter 2. It can be seen that for all 
the models, both with and without secondary heat exchanger, reactor stability is indicated 
at the normal fullpower steady state operation. This is also the case for other reactors 
studied. This is indicated in Table 6 using six groups of delayed neutrons. 


Table 5: Reactor Period for operation at full power for the Reference LMFBR 
(SUPER PHENIX) 


Reactor Model 

Negative Reactor Period (s) 

Without HE 

With HE 

No delayed neutrons 

Effective life-time model 
Prompt-jump approximation 

One group of delayed neutrons 
Six groups of delayed neutrons 

1.31 X 10^ 

1.02 xl0‘ 

2.57 X lO' 

2.57 X 10‘ 

9.15x 10^ 

1.22 X 10^ 

0.8 X 10* 

0.45x10* 

2.42 X 10* 

9.1 lx 10* 
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Table 6: Reactor periods at full power operation for the LMFBRs shown in Table 1 
(using six groups of delayed neutrons) 



Negative Reactor Periods (s) 

jcvcdC/ior 

Without HE 

With HE 

SUPER PHENIX 
BELOYARSK-3 
BELOYARSK-4 

DFR 

EBR-2 

JOYO 

KARLSRUHE 

MANGYSHLASKI 

MONJU 

PHENIX 

FBTR 

9.11 X 10‘ 

9.10x 10‘ 
9.10x10’ 

9.23 X lO’ 
9.83x10’ 

9.04x lO’ 

9.02 X lO’ 
9.14x10’ 

9.06 X lO’ 

9.06x lO’ 

9.15X lO’ 

9.15X lO’ 

9.12X lO’ 

9.12X lO’ 

9.32 X lO’ 

10.5 X lO’ 
9.08x10’ 

9.05 X lO’ 

9.20 X lO’ 
9.10x10’ 
9.10x10’ 

9.12X lO’ 


4.2 Critical Parameter Values 

For the purpose of bifurcation studies, we have treated the coolant temperature coefficient 
of reactivity, , as the variable parameter, keeping all other parameters fixed. (A 
simultaneous variation of all parameters with a view to find the worst combination is 
presented later.) The critical values a/ for the static and dynanaic bifurcations for the 
reference reactor SUPER PEENIX core are presented in Table 7 using 

• No delayed neutrons 
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• Effective lifetime model 


• Prompt-jump approximation 

• One group of delayed neutrons 

• Six groups of delayed neutrons 


It can been seen that the critical values for Hopf bifurcation differ considerably on the 
model employed. In fact the, prompt-jump approximation does not yield any critical 
values for the Hopf bifurcation. For this reason, unless otherwise indicated, we have used 
one-group and six-group models only. 

The effect of secondary heat exchanger on the critical values, a^, is brought out in 
Table 8 and 9, for the transcritical and Hopf bifurcation, respectively, using one and six 
groups of delayed neutrons. While the effect of the number of groups of delayed neutrons 
is relatively small, the effect of the secondary heat exchanger on the critical value is 
considerable. 

Table 10 shows the critical values for other reactors using six groups of delayed 
neutrons and including the secondary heat exchanger in the model. It can be seen from 
Tables 7 to 10, that the static (transcritical) bifurcation occurs for positive value of the 
coolant (sodium) temperature coefficient of reactivity (aj, while the dynamic (Hopf) 
bifurcation occurs for negative values of the same. For most LMFBRs, the acmal values 
of ttj is likely to be negative, and of the order of 10'^ K'\ While the order of magnitude of 
ttc* (transcritical bifurcation) is close to it, it is opposite in sign. On the other hand, while 
the sign of (x^ for Hopf bifurcation is negative, its magnitude is extremely high. Thus 
Hopf bifurcation in any of the reactors studied is completely ruled out. A static 
bifurcation may occur if the coolant (sodium) temperature coefficient of reactivity 
becomes positive. 
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Table 7; Critical values of the coolant temperature coefficient of reactivity (a J for 
Static and Dynamic bifurcations for the Reference LMFBR core (SUPER PHENIX) 


Reactor Model 

Critical value , a g K* 


Static 

Dynamic 

No delayed neutrons 

4.88x10'^ 

-6.14x10'^ 

Effective life time model 

4.88 X 10’^ 

-3.99x10^ 

Prompt jump approximation 

4.88 X 10'^ 

-3.29 

One group of delayed neutrons 

4.88x10^ 

Six groups of delayed neutrons 

4.88x10'^ 

-3.13 


Table 8: Effect of Secondary Heat Exchanger on the critical values for Static Bifurcation 
for the Reference LMFBR (SUPER PHENIX) 


Reactor Model 

Critical Value, * (K'^) 

Core 

Core + HE 

One group of delayed neutrons 

Six groups of delayed neutrons 

4.88 X 10'^ 

4.88 X 10’^ 

3.09 xlO'"' 

3.09 X 10'' 



Table 9. Effect of Secondary Heat Exchanger on the critical values for Dynamic 
Bifurcation for the Reference LMFBR (SUPER PHENIX) 



Critical Value, a/ (K'‘) 

Reactor Model 



Core 

Core + HE 

One group of delayed neutrons 

-3.29 

-4.51 

Six groups of delayed neutrons 

-3.13 

-4.36 


Table 10: Critical values of coolant temperature coefficient of reactivity (a^) for Static 
and Dynamic Bifurcations for the LMFBRs listed in Table 1. 

(Core + Secondary Heat Exchanger using six groups of delayed neutrons) 


Reactor 

Critical values of K*’ 

*( Static) 

a ^ *( Dynamic) 

SUPER PHENIX 

3.09 X 10’^ 

-4.36 

BELOYARSK-3 

3.01 X 10'^ 

-4.44 

BELOYARSK-4 

3.01 X 10'^ 

-3.13 

DFR 

2.02 X 10'^ 

-1.61 

EBR-2 

1.09x10'^ 

-0.836 

JOYO 

4.37 X 10'^ 

-2.78 

KARLSRUHE 

3.87 X 10'^ 

-12.0 

MANGYSHLASKI 

2.78 X 10’® 

-1.07 

MONJU 

4.02 X 10'^ 

-7.44 

PHENIX 

3.54 X 10'^ 

-6.61 

FBTR 

3.01 X 10'^ 

-2.46 
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4.3 Simultaneous variation of Parameters 


As we have seen above, in each of the reactors studied, bifiircation occurs for possibly 
unrealistic values if parameters other than are considered as constant . In this section, 
we vary other dimensionless parameters in certain ranges, with a view to find those sets 
of parameters for which bifurcations may occur for values of a.^ more likely to occur in 
practice. 

For this purpose the values attained by the dimensionless parameters for the SUPER- 
PHENIX reactor are treated as reference values and denoted by the subscript R. These are 
reproduced in Table 1 1 along with the ranges explored for each parameter. Ranges are 
selected with the help of Table 4, with a view to allow for the worst-case combination. 
The actual ranges tried, as given in Table 11, are broader than those strictly suggested in 
Table 4. 

The worst case combinations of parameters are found using random trials. The results 
are shown in Table 12, for both static and Hopf bifurcations, together with the 
corresponding values of a^* ( with and without secondary heat exchanger ) using the 
expressions for the dimensionless parameters a^,b,c,0,p,r,Q given in chapter 2, the 
following physical interpretation may be given: both types of bifurcation are favored by 

• A small magnitude of af 

• A larger value of A and a small value of p 

• A larger value of fuel to coolant heat capacity ratio 

• Operating at larger power P 

• A higher operating fuel temperature 

• A lower heat capacity on the secondary side of the heat exchanger 

Since most of the physical parameters enter into more than one dimensioless variables the 
above observations emphasize only the dominant effects. 



Table 1 1 : Ranges of Parameter Ratios allowed for the worst case LMFBR 



:e Value^ 

Parameter Ratio 

Range 

Refereni 

Lower Limit 

Upper Limit 


-1.045 


2.5 

0.25 

ba 

2x10-^ 

b/ba 

0.1 

10.0 

Cr 

0.652 

C/Cr 

0.4 

3.0 

0R 

0.081 

O/Or 

0.5 

3.0 

Pr 

0.000019 

P^Pr 

0.01 

10.0 

Tr 

0.000029 

r/p 

0.2 

5.0 

©R 

0.667 

©/0R 

0.75 1 

1.5 


t Reference Values are for SUPER PHENIX as given in Table 2 and 4 and are here 
denoted by a subscript R 


Table 12: Values of the Parameter Ratios for the worst case LMFBR for Static 
and Dynamic Bifurcations within the limits given in Table 1 1 



Worst case for 

Parameter Ratios 

Static Bifurcation 

Dynamic Bifurcation 

af / afR 

0.25 

0.25 

b/bR 

10.0 

10.0 

C / Cr 

3.0 

3.0 

O/Or 

0.5 

0.5 

p/PR 

10.0 

10.0 

t/Fr 

5.0 

5.0 

©/©R 

0.75 

1.5 

ac’( Without HE) 

+1.22xlO‘^K'‘ 

-2.6392 X 10"’ K‘‘ 

a; (With HE) 

+7.42xl0‘^K'‘ 

-3.8924 K‘‘ 
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It is, however, to be noted again, that even in the worst case situation, bifurcation values 
are impractical, except when the sodium temperature coefficient of reactivity mav 
become positive. 


4.4 Trajectories, Bifurcation Diagrams and 
Phase Portraits 


As we have pointed above, an LMFBR dynamical system is unlikely to experience a 
bifurcation with the actual value of the parameters. In this section, we present 
numerically simulated trajectories/orbits and the bifurcation diagram for the worst case. 
Furthermore, as has been mentioned in Section 3.3, the static bifurcation is possible only 
when higher-order feedback effects are present. In the figures that follow, these effects 
are included for the transcritical bifurcation. In the case of Hopf bifurcation, we have 
ignored the higher-order feedback effects. 


Figures 5 to 7 shows the time response for the dimensionless power p, 
dimensionless fuel temperatute , and dimensionless coolant temperature 3^ 

subsequent to transcritical bifurcation. The parameter values for which these figures are 
plotted are shown on the figures. It can be seen that while the behavior of the linear 
system, or the system with linear feedback, is divergent, as expected, the behavior 
becomes bounded, as soon as any of the higher-order feedback effect is included. The 
presence of the new equilibrium point is indicated by the asymptotic approach for the 
same in Figures 5 to 7. The actual values at the new equilibrium point depend on how far 
the parameter ac is from the bifurcation point. This is shown in Figure 8 to 10, 
where 


a —a 

Epsilons 



The exchange of stability between the operating point and the new fixed point at the 
bifurcation is also clearly visible in these figures. The solid curves represent the stability 
while the dashed curves represent unstable condition. 

It should however be mentioned that the new equilibrium point is only locally 
stable( in accordance with the local bifurcation theory ). Thus all disturbed situations 
from the original operating point need not be in the basin of attraction of the new fixed 
point. The behavior shown in Figures 5 to 7 is valid only for those disturbances which are 
in the basin of attraction of the new equilibrium point. For other disturbances, the reactor 
behavior will remain divergent. 

Figures 11 to 15 show the behavior subsequent to Hopf bifurcation in an LMFBR 
dynamical system using one and six groups of delayed neutrons, in the worst case. As has 
already been mentioned, the value of the parameter for which this behavior occurs 
appears to be extremely impractical. These figures, however, confirm the predictions of 
the theory described in Chapter 3. It may also be noted, that even this unlikely situation 
while the neutron flux may show large peaks, their effect on fuel and coolant temperature 
is limited ( of the order of 20% to 30%). 

The stability of the limit cycles represented in Figures 11 to 15 is further confirmed 
by the calculation of stability parameter, a(0), (3.2) and the characteristic multipliers 
shown in Table 13 . We mention here that as a(0)>0 , the situation similar to that shown 
in Figure 4c occurs in these models. In view of the unphysical nature of , we will not 
consider the unstable limit cycles in the vicinity of the bifurcation point. These and their 
associated effects have been investigated by Manmohan (1996) in great detail. 
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4.5 Critical Values for the Control 
Parameters 


As given in section 2.9, there are four control parameters Kp, tp, K„ and t,. Table 14 
gives the critical values oc^ , for different values of Kp and tp as compared to the situation 
without control ( Kp— 0, tp = 0 ). It can be seen that the bifurcation becomes all the more 
unlikely for all values of Kp < 0, the only situation that can occur in practice. 

Table 15 shows the critical values of K^ for different values of t^. It can be seen that 
Kc* is of the order of (- 10^ to -10^ ). This, for the values of T<,o, Tjo, and p given in Tables 
1 and 2, a control reactivity feedback of the order of -100 to -1000 $ per degree rise in 
sodium temperature. A practical value for the same is likely to be of the order of -10’^ to 
-10’ $ per degree rise in the coolant temperature. Thus it can be concluded that the 
control system is unlikely to introduce any kind of bifurcation in the system. 

Table 13: The stability parameter a and the characteristic (Floquet) multiplier for the 
LMFBR core limit cycle( worst case) 


Parameter 

One Group 

Six Groups 

Stability Parameter , a(0) 

Floquet Multipliers (s = 0.005) 

6.421x10"^ 

0.999 

0.235x10'' 

0.512x10'" 

6.426x10'" 
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Table 14 : Critical values of the coolant temperature coefficient of reactivity (Uc) for 
different values of Kp and tp for the reference reactor (SUPER PHENIX) Core + HE using 
six groups of delayed neutrons 



tp(s) 

Critical values, 






Static Bifurcation 

Dynamic Biftircation 

0.0 


3.09x10'^ 

-4.36 

-0.1 

0.0 

8.31x10-^ 

-3.7x10^ 


0.01 

2.96x10'^ 

-4.126 


0.1 

0.288 

-3.64 

-1.0 

0.0 

8.03x10'^ 

.3.48x10^ 


0.01 

2.96x10-' 

-3.99x10' 

-10.0 

0.0 

8.00x10'^ 



Table 15: Critical values of the control parameter IC,. (Kp=0, tp=0) for Hopf biftixcation for 
Reference LMFBR (SUPER PHENIX). 


tc(s) 

Kc* 

(dimensionless) 

0 

-2.18x10' 

LOxlO-’ 

-3.35x10' 

2.0x10’’ 

-7.26x10' 

2.5x10-’ 

-1.74x10' 


central library 

« »■ T., kanw 


4.6 Conclusions and Observations 


Based on the analytical and numerical results presented in this work, the following 

conclusions can be drawn : 

• Models based on effective lifetime and prompt-jump approximations do not yield 
correct critical values. While these models may be quite useful in linear analysis, 
these should not be used for bifurcation studies. 

• In lumped parameter dynamical systems including the core and the secondary heat 
exchanger, based on one and six groups of delayed neutrons, Hopf bifurcation is ruled 
out for physically realistic values of the parameters. This remains true for all 
(negative) gains in the proportional controllers, and any time constants associated 
with the differential controllers. 

• Static bifurcation can occur only if the higher order feedback coefficients are present, 
of the three types of static bifurcation, saddle-node, pitch-fork and transcritical 
bifurcations it is shown in the present work that in the models considered only 
transcritical bifurcation can occur. However, this can occur only for positive values of 
coolant temperature coeflEicient of reactivity (a^ ) of the order of 10’^ K‘\ For most 
LMFBRs, the values of ac, while of the same order of magnitude, is likely to be 
negative. 

• It is noted in the present work that subsequent to transcritical bifurcation, the new 
equilibrium point is only locally stable and may not be an attractor for all 
disturbances from the normal operating conditions. Thus there remain disturbances of 
the operating condition from which the reactor will show unstable divergent behavior 
in spite of the presence of higher order feedbacks 

Some of the limitations of the present work and a few observations regarding future 

work in this area are given below: 
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. We have used lumped parameter models in the present work. Inclusion of the space 
variables in the models can modify / change the above conclusion. This should be 
included in future studies. 

• In the present work, the control reactivity, whether based on change in reactor power 
or based on change in the coolant temperature, is assumed to act instantly. In future 
studies, appropriate governing equations for the control reactivity insertion should be 
included. 

• We have ignored the transport times in the coolant pipes external to the reactor, 
leading the coolant from the reactor to the secondary heat exchanger. If the control 
signal based on coolant temperature is derived from measurement of the same across 
the secondary heat exchanger, the above mentioned transport times can play an 
important role. This should be included in future studies, together with any heat loss 
in the coolant pipes external to the reactor. The later effect will also modify the heat 
exchanger model. We may mention that the inclusion of the transport delay in 
nonlinear models necessitates a consideration of bifurcation theory of delay 
differential equation which has been outside the scope of the present work. 

• In the present work, we have studied the local behavior in the vicinity of the 
bifurcation points. A global analysis may be attempted in future studies. 

• Nonlinearities also effect the transient behavior of the reactor under normal operating 
parameter values. Their effect on the transients may have an important bearing on the 
optimal design of the controllers / control law. This can be a fruitful exercise to 
undertake in future. 
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Dimensionless Power 
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Figure 5: Time variation of dimensionless power p subsequent 
to transcriticai bifurcation in an LMFBR core (worst case) . 
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Figure 6 : Time variation of dimensionless coolant temperature 
3, subsequent to transcritical bifurcation in an LMFBR core 

(worst case). 
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Figure 7: Effect of secondary heat exchanger and groups of 
delayed neutrons on time variation of dimensionless fuel 
temperature 3 , subsequent to transcritical temperature in a 

LMFBR core. 
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Figure 9 : Bifurcation diagram (with and without secondary 
heat exchanger ) for transcritical bifurcation in the worst case 

LMFBR 


/ 

[cCf etc 
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Figure 10 : Bifurcation diagram (with and without secondary 
heat exchanger ) for transcritical bifurcation in the worst case 
LWIFBR 
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Time(s) 


Figure 11 : Behavior of Dimensioniess Power (neutron flux ) 
versus time in LMFBR core (worst case) subsequent to Hopf 
bifurcation using one group of deiayed neutrons ( s = o.oi ) 
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Dimensionless Power 


Figure 12 : Phase portrait of limit cycles in LMFBR models on 
plane (worst case) subsequent ^ to Hopf bifurcation 

(^== 0 . 01 ). 
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Dimensionless Coolant Temperature 



Figure 13 : Phase portrait of limit cycles in LMFBR models on 
3/ - plane (worst case) subsequent to Hopf bifurcation 

(^ = 0 . 01 ) 
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Time(s) 


Figure 14 ; Behavior of Dimensionless Power (neutron flux ) 
versus time in LMFBR core (worst case) subsequent to Hopf 
bifurcation using six groups of delayed neutrons (s = o.oi ) 
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Figure 15 : Effect of secondary heat exchanger on the limit 
cycles In the worst case LMFBR subsequent to Hopf bifurcation 

(6“«0.0l) 
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Appendix A 


Sample Calculation of Cf and Cc for 
Reference LMFBR SUPER PHENIX 

For this reference reactor SUPER PHENIX, we have, (Power Reactors (1983)) 

Estimation of Cf 

Active core height, H = l.Om 

Active core diameter D = 3.66m 

Fuel pellet diameter, d =7.14xl0"^m 

Clad thickness, t =.57xl0"^m 

Fuel inventory, mf = 31.5 x 10^ kg 

Fuel density (U/Pu metal), pf = 19000 kg m'^ 

Specific heat (U/Pu metal), Cf = 0.17 KJ kg'* K'* 

Cf=mfCf= 5.355 MJK'* 

Estimation of Cc 

Core volume, V = (7t/4)D^H = 10.62 m^ 

2 

Fuel volume , Vfoei = nif/ Pf = 1.6 m 

69 


^structure ^ ^fiiel C^clad ^ ^fiiel) X 1.2 — (t(d+t) / ) xl.2 —0.414 

(20 /o additional structure volume over and above the clad volume allowed in the core.) 
Vstructure ' V = 0.16 x0.414 = 0.066 

(VfUel + V3„e) / V = 0.226 

Vcoolant/V = 0.774, 

VcooIant=8.21m^ 

Pcooiant (Na) = 800 Kgm'^ 
nicooiant= PcVc = 800 x8.21 = 6576 Kg 
Cc = mcC, = 7.91MJK'‘ 

Other Reactors 

The core/fuel data for other reactors is given in Table A. For calculations in the present 
work , densities and specific heats are taken as follows 

p(U02/Pu02)= 11000 kgm’^ 
c (UO 2 / PUO 2 ) = 0.32 KJ kg'* K‘‘ 
p(UC/PuC) = 13630 kg m'^ 
c (UO 2 / PUO 2 ) = 0.25 KJ kg'*K'’ 
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Table A : Core fuel data for the representative LMFBRs 


Reactor 

Fuel 

Inventory 

(t) 

Core 

Height 

(m) 

Core 

Diameter 

(m) 

Fuel Rod 

Diameter 

(m) 

Clad 

thickness 

(m) 

SUPER PHENIX 

31.5 

1.0 

3.66 

X 

O 

0.57. 10'^ 

BELOYARSK-3 

8.5 

0.75 

2.05 

5.9.10'^ 

0.4. 10'^ 

BELOYARSK-4 

8.5 

0.75 

2.05 

5.9. 10'^ 

0.4.10"^ 

DPR 

0.34 

0.53 

0.53 

4.95.10'^ 

0.5.10'^ 

EBR-2 

0.46 

0.36 

0.65 

4.4. 10’^ 

0.3x10’^ 

JOYO 

0.25 

0.55 

0.73 

4.6. 10'^ 

0.35.10'^ 

KARLSRUHE 

0.728 

0.6 

0.82 

6.4. 10'^ 

0.5x10'^ 

MANGYSHLASKI 

1.22 

1.0 

1.58 

5.9.10'^ 

0.4x10'^ 

MONJU 

5.9 

0.93 

1.79 

o 

X 

X 

O 

PHENIX 

4.3 

0.85 

1.39 

5.1.10’^ 

0.45.10'^ 

FBTR 

0.19 

0.32 

0.43 

3.65.10'^ 

0.37x 10'^ 
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